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ABSTRACT 

Massive neutrinos make up a fraction of the dark matter, but due to their large 
thermal velocities, cluster significantly less than cold dark matter (CDM) on small 
scales. An accurate theoretical modelling of their effect during the non-linear regime 
of structure formation is required in order to properly analyse current and upcoming 
high-precision large-scale structure data, and constrain the neutrino mass. Taking ad- 
vantage of the fact that massive neutrinos remain linearly clustered up to late times, 
this paper treats the linear growth of neutrino overdensities in a non-linear CDM 
background. The evolution of the CDM component is obtained via iV-body computa- 
tions. The smooth neutrino component is evaluated from that background by solving 
the Boltzmann equation linearised with respect to the neutrino overdensity. CDM 
and neutrinos are simultaneously evolved in time, consistently accounting for their 
mutual gravitational influence. This method avoids the issue of shot-noise inherent 
to particle-based neutrino simulations, and, in contrast with standard Fourier-space 
methods, properly accounts for the non-linear potential wells in which the neutri- 
nos evolve. Inside the most massive late-time clusters, where the escape velocity is 
larger than the neutrino thermal velocity, neutrinos can clump non-linearly, causing 
the method to formally break down. It is shown that this does not affect the total 
matter power spectrum, which can be very accurately computed on all relevant scales 
up to the present time. 

Key words: neutrinos - cosmology: large-scale structure of Universe - cosmology: 
dark matter 



1 INTRODUCTION 

The determination of neutrino masses lies on the bound- 
ary between particle physics, astrophysics and cosmology. 
Atmospheric and solar neutrino oscillations have allowed 
the measurement of two mass-squared differences AmJ 2 = 
ml-ml = 7.54±o.39 x 1CT 5 eV 2 and | Am 2 3 | = |m| - (m 2 + 
m!)/2| = 2.43 1" 2 , x 10 ~ 3 ey2 (central values and 2-a er- 
ror bars from |Fogli et al.||2012| . These measurements do 
not allow us to pin down the individual neutrino masses, 
nor their hierarchy, since only the absolute value of Am 2 3 is 
measured - note, however, that next-generation oscillation 
experiments such as NO\AQaim to distinguish the neutrino 
hierarchy through second-order mixing effects. Current /3- 
decay experiments imply an upper limit on the individual 
neutrino mass m„ ^ 2 eV ( Kraus et al.|2 005), and future ex- 
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periments like KATRIN should lower this limit by an order 
of magnitude ( |Eitel|2005| ). 

Cosmological probes are sensitive to the absolute mass 



of neutrino species (see, for example, Lesgourgues & Pastor 
2006 and Wong|20"TT for recent reviews), and so far provide 
the tightest limits on the total neutrino mass M v = ^ rn v = 
mi + m2 + rn-i, at the level of a few tenths of electron- 
volts (eV). Increasingly sensitive future experiments aim 
to reach levels of precision at which a total neutrino mass 
~ 0.1 eV or lower would be detectable (we refer the 
reader to Abazajian et al. 2011 for a recent compilation of 
current constraints and forecasts). Ultimately, the neutrino 



hierarchy may even be within reach ( |Takada et al. 
Jimenez et al.|2010[ ). 



2006 



NOVA proposal: http://www-nova.fnal.gov/ 



Massive neutrinos affect the growth of perturbations 
through two effects. First, they change the background evo- 
lution. For a total neutrino mass less than a few eV, neutri- 
nos are still relativistic at the time of matter-radiation de- 
coupling but non-relativistic today. At fixed total density of 
non-relativistic matter today, the epoch of matter-radiation 
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equality is moved to later times as the neutrino mass is in- 
creased. If one holds the density of baryons and CDM fixed 
instead, the angular-diameter distance and time of matter 
- dark energy equality are changed for different neutrino 
masses. In each case the CMB anisotropy power spectrum 
is affected ( |Lesgourgues fc Pa stor 2006; Wong 20TT|, which 
has allowed CMB observations to set what is perhaps the 
most robust constraint on the neutrino mass, ^ ra„ < 1.3 
eV ( |Komatsu et al.||2011[ ). Upcoming data from the Planck 
satellite is expected to lower this bound by a factor of two 



few percent accuracy up to k ~ 1 h Mpc at high red- 



I Planck Collaboration 2005 Abazajian et al. 20111 



The second effect of massive neutrinos is to slow down 
the growth of structure on scales smaller than their free- 
streaming length. This effect has long been understood in 



the linear regime (Bond et al. 1980 Bond & Szalay 1983 Ma 



fc Bcrtschingcr 1995 . Lcsgourgucs & Pastor 2006 , and refer- 
ences therein) , and is implemented in all modern Boltzmann 
codes (see for example Lewis fe Challinor|2002 |Lesgourgues 



fc Tram|2011 1. For a total neutrino mass of a few tenths of 
eV, this effect becomes manifest at scales £ 100ft -1 Mpc 
and gets stronger in the non-linear regime at scales of a few 
tenths to a few tens of Mpc ( Lesgourgues fe Pastor|20"06l ). 

Combined with CMB anisotropy measurements, various 
probes of the matter distribution have already constrained 
the total neutrino mass to ^ra u 0.2 — 0.3 eV. These in- 



clude galaxy surveys such as SDSS (do Putter et al. 2012 Xia 
et al. ||2012h , cosmic shear surveys such as CFHT-LS (Ichiki 



et al. 



|2009|), Lyman- q forest measurements (Seljak et al 
Viel et al.|2010 1, and galaxy cluster surveys ( Vikhlinin 



20061 

et al. 20091. As the survey samples grow larger, the sensi- 
tivity to the neutrino mass is expected to reach ^ m„ ~ 0.1 



eV ( Abazajian et al7||2011 1 or even lower (see |Takada et al 
2006 for a detailed Fisher-matrix forecast of the constraining 
power of future high-redshift surveys). In addition, future 
measurements of CMB lensing ( Hall fc Challinor|2012 I and 
21 cm surveys (Ma o et al.|2008 1 have the potential to reach 
yet lower bounds. Massive neutrinos also affect the growth 
of halos ( |Brandbyge et al.|2010[ |Marulli et al.|2011[ ) . 

Each one of the aforementioned methods has its own ad- 
vantages and limitations. Systematic errors may arise from 
using biased tracers of the underlying density field, from 
difficulties in modelling complex baryonic physics, or from 
foreground contamination. In addition, this level of precision 
requires an accurate modelling of corrections from non-linear 
structure growth. Here we focus on the latter problem. Cos- 
mological neutrinos only interact gravitationally and, be- 
cause of their large thermal velocities, effectively constitute 
a hot dark matter (HDM) component. While it poses signif- 
icant practical difficulties, the non-linear clustering of mas- 
sive collisionless particles is now a relatively well-understood 
problem. The growth of pure cold dark matter (CDM) struc- 
ture can be computed through collisionless iV-body simula- 
tions, whose accuracy is limited primarily by available com- 
putational resources, which have been steadily increasing 



over time (see e.g. Fig. 1 of Dehnen fc Read|2011[ ). By con- 
trast, the effect of neutrinos on the non-linear growth of 
matter perturbations has only recently been investigated in 
depth. This delicate problem has been approached from two 
angles. On the one hand, extensions to perturbation the- 
ory allow for semianalytic calculations of the CDM power 
spectrum in the weakly non-linear regime. Combined with 
a linear treatment of neutrinos, such methods can reach a 



shifts ( z £ 2) ||Saito et al.||2008| |2009| |Lesgourgues et al 



[2009}|Wong||2008| ). |Shoji fc Komatsu| ( |2009| ) go beyond lin- 
ear theory for neutrinos and solve for both the CDM and 
neutrino overdensities to third order in perturbation theory, 
making however the simplifying assumption that neutrinos 
can be described by an ideal fluid with a constant Jeans 
scale. Hannestad et al. (20121 incorporated neutrinos into 



N-body simulations using a similar fluid approach. 

On the other hand, the "exact" solution to the prob- 
lem of CDM+neutrino clustering can in principle be ob- 
tained from iV-body simulations where both CDM and neu- 
trinos are simulated as particles. This approach was recently 
taken byjBrandbyge et al.| ( |2008[ ), |Viel et al.| ( |2010[ | and |Bird| 
et al.| ( |2012[ ) (see also references therein for earlier works). 
The main limiting factor of such an approach is numeri- 
cal: since neutrinos do not significantly cluster below their 
free-streaming scale, their power spectrum on small scales is 
dominated by shot noise for any reasonable number of sim- 
ulated particle^] increasingly so for lower neutrino masses 
and at earlier times. A vast amount of computational power 
and memory is therefore wasted to extract a small amount 
of information on the actual neutrino clustering. In order to 
bypass the shot-noise issue, |Brandbyge fc Hannestad] ( |2009[ ) 
included neutrino perturbations in Fourier space, while still 
treating the CDM component as particles in iV-body simula- 
tions. The density field of their neutrino component was only 
computed using linear theory, however, and did not account 
consistently for non-linear CDM perturbations sourcing the 
gravitational potential in which neutrinos evolve ( |Bird et al.| 
20121. Recently, Brandbyge & Hannestad (20101 suggested 



a hybrid method, where neutrinos are initially treated us- 
ing the Fourier method with linear theory, and converted to 
particles as time progresses. One of the remaining sources of 
errors of this implementation is that the Fourier-space part 
does not account for the non-linear growth of gravitational 
potentials in the neutrino evolution. 

The aim of the present work is to close this gap and pro- 
vide an efficient method for accurately computing the effect 
of massive neutrinos on the non-linear matter power spec- 
trum, with a minimal modification to well-tested pure CDM 
iV-body simulation methods. To do so, we follow the CDM 



with the Tree-PM iV-body code Gadget-3 ( Springel||2005 



Viel et al. 20101, and compute analytically the linearised 



neutrino overdensity sourced by the full gravitational po- 
tential. The CDM and neutrinos are evolved simultaneously 
and self-consistently, accounting for their mutual gravita- 
tional influence. Our method is therefore semi-linear, in the 
sense that we account exactly for the non-linear growth of 
CDM overdensities and gravitational potentials, but effec- 
tively use a linear transfer function to obtain the neutrino 
overdensity from the gravitational potential. For neutrinos 
of the mass allowed by current data, our implementation is 
sufficient on its own for a complete description of the matter 
power spectrum up to the present time, and, for z > 1, the 



2 The shot noise issue arises for neutrinos for two reasons. First, 
their large random velocities effectively randomly distribute them 
over the box, which leads to a Poisson spectrum P{k) = 1/n at 
all scales. Second, their intrinsic clustering is very small on small 
scales. At large enough redshifts and for small enough scales, the 
true power may be completely swamped by the Poisson noise. 
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neutrino power spectrum as well. It only fails to resolve neu- 
trino clustering in the most massive galaxy clusters, in which 
neutrinos do cluster non- linearly (Ringwald & Wong 2004[ ); 
however, it should be easy to perform zoomed simulations of 
these clusters using our method as a base. An ancillary ad- 
vantage of our implementation is that it can easily be added 
into an TV-body code as a small "patch" and does not re- 
quire running an independent Boltzmann code in parallel. 
Furthermore, our code can easily include the exact contri- 
bution of neutrinos to the background expansion, and the 
effect of the neutrino hierarchy, which are problematic or 
expensive to include in particle implementations. It can be 
applied to regimes where shot noise severely limits the reli- 
ability of the particle method, such as low neutrino masses, 
the 21cm forest, or large- volume galaxy mock catalogue sim- 
ulations. Perhaps the most important property of our code 
is that it allows massive neutrinos to be included into sim- 
ulations with negligible cost in CPU and memory, and with 
no loss of accuracy. Thus it can be useful for interpreting 
the results of essentially any probe of large scale structure, 
including CMB lensing, the Lyman-a forest, weak lensing 



experiments or galaxy surveys ( Takada et al.||2006 Cooray 
19^91 |Kaplinghat et al.||2003| |Wang et al.||2005| |Vallinotto[ 
et al'||2009| |Gratton et al.||2008| |Ichiki et al.||2009| | Schlegei| 



et al.|(2009[ Carbone et al. 2011 |Beaulieu et al.|2010[ >. 

This paper is organised as follows. In Section [2] we de- 
fine our notation and discuss characteristic scales and the 
regime of validity of our main approximation. We derive our 
main equations in Section [3] The practical interfacing with 
our TV-body code is detailed in Section]!] We describe our re- 
sults and compare them against other methods in Section [5] 
We discuss potential applications of our method in Section 
[6] and conclude in Section [7] 



2.2 Relativistic to nonrelativistic transition 

We denote by p the physical momentum of a neutrino and 
fo (p) the unperturbed neutrino phase-space density, giving 
the number of neutrinos and antineutrinos per unit physical 
phase-space in the absence of gravitational clustering. Neu- 
trinos decouple from the baryon plasma at a temperature 
of about 1 MeV (ILesgourgues & Pastor 2006 1 . Given that 



their masses are known to be less than 1 eV, this means 
that they were ultrarelativistic at decoupling. This implies 
that their unperturbed phase-space density is the redshifted 
relativistic Fermi-Dirac distribution (since p oc a -1 in the 
absence of gravitational potentials), 



/o(p) = 



1 



/i 3 exp(p/T„) + l' 



(3) 



with T v (z) — (1 + z)T„,o. A neutrino of mass m„ becomes 
non-relativistic when its momentum falls below its mass, 
p & m v . This occurs at a redshift z nl such that 



m 



^ ft .-,<).-,- 



(4) 

0.1 eVp w 

The Fermi-Dirac distribution |3| is such that 50% of neu- 
trinos have a momentum p < 2.84 T v and 90% have a mo- 
mentum p < 5.47 T„. As a consequence, 50% of neutrinos 
are nonrelativistic for z ^ 200 (m„/0.1 eV) and 90% have 
become non-relativistic by z ^ 100 (m„/0.1 eV). For the 
lowest masses considered (m v « 0.05 eV), relativistic cor- 
rections to neutrino clustering may be of order unity at high 
redshifts (we start our simulations at z = 49); however, neu- 
trinos are basically unclustered on all scales when they are 
quasi-relativistic, and the exact value of their very small 
inhomogeneities is then irrelevant for CDM clustering. We 
therefore treat neutrinos in the non-relativistic limit, i.e. as- 
sume p <C m„ in all our derivations^] 



2 NOTATION, CHARACTERISTIC 

TIMESCALES AND LENGTHSCALES 

2.1 Notation 

Throughout this paper we use r for the conformal time 
(dr = dt/a, where t is the physical time and a = (l + z)~ 1 is 
the scale factor) and work in units where c = G = ks = 1. 
Lengthscales and wavenumbers are comoving unless other- 
wise stated. Overdots denote differentiation with respect 
to conformal time and gradients are with respect to co- 
moving lengths. The Hubble rate today is Ho = 100 h 
km/s/Mpc. When referring to an individual neutrino mass, 
we use the lower-case m v . When referring to the sum of 
neutrino masses, we use the upper-case M v = 5]m„. The 
temperature of the (unperturbed) neutrino background to- 
day is 



T„, = 1.95 K = 1.68 x 10" 4 eV. 



(1) 



Neutrinos being non-relativistic today, they contribute a 
fraction /„ of the total nonrelativistic matter abundance, 
given by 



1 



n M h 2 93.14 eV 



0.022 



0.147 £m„ 
n M h 2 0.3eV 



(2) 



2.3 Free-streaming scale 

Neutrinos can free-stream across a comoving lengthscale A 
if the time it takes them to cross A is much less than the 



Hubble time H 



i.e. if 



aX 



(5) 



where v v (z) ~ T v {z)/m l , is the characteristic neutrino ve- 
locity. This equation defines a characteristic comoving free- 
streaming scale kf s (z) ~ aH/v v (z). A more detailed analysis 
in Section[3]will show that it is convenient to define the free- 
streaming scale as 



fcf s (a) 



-Q. M (a)v- 



1/2 



(6) 



where v~ 2 is the mean inverse velocity squared, 



21n(2) 



m v a 



810 km/s (1 + Z )°^-?X ] 5 ( 7 ) 



3 Note that consistently accounting for relativistic corrections in 
the evolution of neutrinos would also require accounting for CMB 
inhomogeneities sourcing the gravitational potentials; these are 
always (rightfully so) neglected in TV-body simulations. 
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Figure 1. Non- linear scale (solid) and free-streaming scales for 
various neutrino masses, as a function of redshift. 



and Am (a) is the relative contribution of the non-relativistic 
components to the total energy density at scale factor a. On 
scales much larger than the free-streaming scale, neutrinos 
behave like CDM. We shall show later on that for fc ^> fcf s , 
the (linearised) neutrino overdensity 5 V relates to the total 
(possibly non-linear) matter overdensity 5m as 



5 v {k > fcf s ) 



Numerically, we obtain 



m " -h Mpc- 1 . 



(8) 



(9) 



We show the non-linear scale fc n i (defined such that the vari- 
ance of CDM overdensity per log-fc interval reaches unity on 
that scale) and free-streaming scale as a function of neutrino 
mass and redshift in Fig. [I] We see that for ^ m v < 0.6 eV, 
we always have fcf s < fc n i, with an increasing difference at 
large redshifts. Moreover, the non-linear matter powe r spec- 
trum typically grows as k 3 P M {k) oc k a , with a <, 2 jSeljak 
|2000[ ) and therefore we expect the power per log interval in 
the neutrino component to be decreasing for fc ^> fc n i > fcf s 
as fc 3 P„(fc) oc fc' 3 , with /3 & -2. 

We therefore expect the neutrino component to be lin- 
early clustered on all scales for masses below the current 
upper limit: for k < fc n i neutrinos cluster at most like the 
CDM, which is itself linearly clumped. For fc > k n \ > fcf s , 
the neutrino power per logarithmic interval is a decreasing 
function of fc. 

This argument motivates our use of linear theory for 
the neutrino component. However, this does not give the full 
physical picture of neutrino clustering: the slowest neutrinos 
may in fact significantly cluster in massive haloes. Before 
moving to the core of our calculation, we first discuss in 
which conditions it may break down. 



2.4 Neutrino capture in massive haloes 

In this section we qualitatively discuss under which condi- 
tions neutrinos may significantly cluster in a potential well 
of characteristic physical scale ro and characteristic depth 
4>o < 0. Let us consider, for simplicity, a spherical top hat 
potential, 4>{r < ro) = <j>o and <j>{r > ro) = 0. We assume 



that ro is constant in time, but allow the depth <j>o to vary 
slowly, on a Hubble timescale, as is the case for the most 
massive halos currently forming. We consider the fate of a 
particle on a purely radial trajectory. If the particle enters 
the potential at time £j, with initial velocity v(t~) — Vi, its 
velocity upon entry becomes 



v(t+) = Jv[+2\MU)\. 



(10) 



Since we have assumed the potential to be flat inside the 
halo, the particle's velocity is conserved until it reaches the 
other end, at time tf. By then the gravitational potential 
has grown a little deeper, and the particle will escape only 
if its velocity is larger than the new escape velocity, i.e. if 

B?+2|0o(«i)|>2|0o(*/)|. (11) 

Provided the crossing time is short compared to the evolu- 
tion timescale of the potential, we can Taylor-expand this 
equation and obtain the escape condition 

B?>2(t / -t t )|^ |. (12) 

Inside the halo, provided 2\<f>o\ 2> vf, the velocity is approx- 
imately y / 2\(f> \ and the crossing time is therefore 

tf -U^^=. (13) 
If we define the timescale for variation of d> as 



At 



(14) 



the escape condition for neutrinos becomes 
r <i(//A^)|-|=, (15) 

where we have purposefully inserted the Hubble parame- 
ter to make the order-unity parameter HAt$ appear. We 
see that for deep potential wells varying on the Hubble 
timescale, the condition for neutrinos to truly free-stream 
and not be captured is more stringent than simply requiring 
their characteristic scale to be smaller than the (physical) 
free-streaming scale rf s ~ v/H. Equation (15 I can also be 
turned into an escape condition for the neutrino momentum. 
For z rs 0, this is 



P_ 



s/H At 4 



2H r vl^oi 



{H At 4 



-1/2 



ro 



0.1 eV \0.5 /j-'Mpc 



1/2 



1/2 



3000 km/s 



1/2 



We have normalised the lengthscale and depth of the grav- 
itational well to values typical for the most massive haloes. 
The outcome of this analysis is that in massive haloes 
varying on a Hubble timescale, neutrinos with momentum 
p S~> T v are typically captured, while those with momen- 
tum p <; T„ can escape. We emphasise that the cutoff value 
at p w T v is a pure numerical coincidence, arising from 
the characteristic sizes and depths of massive haloes, and 
for neutrino masses of order 0.1 eV. Given that only about 
6% of neutrinos have a momentum p < T v , the qualitative 
picture that emerges from this analysis is that a relatively 
small fraction of neutrinos are efficiently bound to massive 
haloes, thereafter strongly clustering, while a majority re- 
main weakly clustered. 
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Let us emphasise that although the crossing time of 
a neutrino in a massive halo is significantly shorter than 
the Hubble time (for the fiducial values of Eq. | |16| , H8t « 
1/40 at 2 = 0), the relevant comparison is not of HSt with 
unity but rather with the small ratio v 2 /<j> (also w 1/40 for 
our adopted fiducial values): even a small relative effect due 
to the change of the potential on a Hubble timescale may 
translate into an order unity effect on the pre-entry velocity 
and prevent the neutrino from escaping the halo. 

The analysis presented here is of course very simplified, 
for we have used an idealised (and unphysical) profile for the 
gravitational potential. Moreover, in reality there is no sharp 
boundary in momentum space between linear behaviour and 
strong clustering. It is however a robust statement to say 
that the characteristic momentum of neutrinos efficiently 
captured in massive haloes is of order of their temperature. 
If the critical momentum were much larger than T„, the bulk 
of neutrinos would be captured and the neutrino overdensity 
in massive clusters would reach values comparable to those 
for the CDM. If it were much smaller than T„, then linear 
theory would be extremely accurate for neutrinos even in 
the vicinity of massive clusters. The actual situation is in- 
termediate between these two regimes, as shown in the work 
of Ringwald & Wong (2004), who used iV-one-body simu- 
lations to track the evolution of neutrinos around massive 
clusters. They found that for m„ — 0.15 eV, linear theory 
underestimates neutrino clustering by a factor of ~ 2 near 
the centre of 10 14 Mq clusters, and by a factor of ~ 4 for 
1O 15 M clusters; for M < 1O 12 M , they found an excellent 
agreement with linear theory for m„ = 0.15 eV. 

With this caveat in mind, we now proceed to the main 
thrust of this paper: linear perturbations of neutrinos in a 
general gravitational potential. 



3 LINEAR EVOLUTION OF HOT DARK 
MATTER PERTURBATIONS IN A 
GENERAL BACKGROUND 

In this section we consider the linear evolution of an en- 
semble of non-relativistic hot dark matter (HDM) particles 
of mass m in a general gravitational potential <f> in an ex- 
panding universe. We derive an integral equation known as 
the Gilbert equation ( Gilbert|1966 1 , which has been used in 
various analytic works on neutrinos ( Bond fc Szalay||1983| 
" Brandenberger et al.||1987| |Singh fc Ma||2003| |Abazajian| 
et al.[ 2005). This equation can easily be generalised to the 
relativistic case jWeinberg|2008[ |Shoji fc Komatsu|2010| [de] 



|Vega fc "Sa nchez 20 12a|b[ ), although it takes on a more com- 
plicated form 



3.1 Vlasov equation 

3.1.1 Derivation for non-relativistic particles in an 
expanding universe 

In principle, the Boltzmann equation in an expanding uni- 
verse should be derived in a consistent relativistic fashion 
(see for example Ma & Bertschingcr 1995). However, for 
non-relativistic particles and in the Newtonian (<j> <C 1) and 
subhorizon (fc 2> aH) limits, which always hold on the scales 
of interest, we can derive it more simply as follows. 



We choose an arbitrary coordinate origin xo and denote 
by r = a(r)(x — Xo) the physical coordinate of a particle, 
where x is its comoving coordinate. We also denote by p the 
physical momentum of the particle, 

dx 



dr 

m— — mhr 

dt 



m — — = rnHv + a 
dr 



(17) 



where the last equality defines the comoving momentum q ; 
am dx/dr. The equation of motion for a test particle is 



-maV x $ 



efa 



(18) 



dp „ T . da 

It = " mVr *' le ' Tr 
where 3> is the total gravitational potential. We now decom- 
pose $ into a piece due to the background density and a 
piece sourced by the density perturbations: 



$ = $o + < 
where 

$o(r) = 



(19) 



1 d 2 a 
~2adt?* 
2tt 



2 V dt 



(p + 3p)r 

and <f> is determined from the Poisson equation 
V 2 T (j> = 4vr5p. 



(20) 



(21) 

Equation ( |18[ ) for the comoving momentum then simplifies 
to 
dq 
dr 



-maV x ( 



(22) 



We now define /(x,q, r) as the HDM phase-space density, 
which is the number of particles per unit of physical phase- 
space (i.e. diVhdm = /(x, q, T)d 3 rd 3 p). In the absence of col- 
lisions, the conservation of phase-space along a particle's 
trajectory gives the collisionless Boltzmann equation (also 
known as the Vlasov equation): 



<1_ 
dr 



tra.j 



0. 



(23) 



Finally, using the definition of q and Eq. ( 22 l for its deriva- 



tive, we find (independently of the chosen origin xo): 

maV x cj) ■ V q / = 0. (24) 



or ma 



3.1.2 Linearization and Fourier transform 

In a homogeneous universe, V x <£ = and V x / = 0. As 
a consequence, the solution to the Vlasov equation must 
satisfy dfo/dr — 0, i.e. be a function of momentum only. 
Finally, isotropy guarantees that /o = fo(q) is a function of 
q = |q| only. 

We now linearise the Vlasov equation about the homo- 
geneous solution: /(x, q, r) = fo(q) +<5/(x,q, r), assuming 
the perturbation Sf <C fo is induced by the presence of the 
gravitational potential cj> (our formal expansion parameter 
is (ma) 2 (f>/q 2 ). The linearised equation becomes (Branden- 
berger et al.|1987| [Ma fc Bertschinger||1995 1 



95 f q _ madfo 

S, 1 • v x o/ p-q ■ Vx<p = 0. 

or ma q dq 



(25) 



We can now Fourier-transform this equation and obtain 
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95 f q-k ma dfo, . 

— h i 0/ = i — -i-(q • k)0, 

or ma g ag 



(26) 



where k is the comoving wavenumber, and 8f(k, q, r) and 
0(k, r) should now be understood as the Fourier transforms 
of Sf and <j>, respectively. 



If we expand the phase-space density at the initial time on 



the basis of Legendre polynomials (see for example Weisstein 
[20124 , 



<5/(k,q,Ti) ='^i L Sfi(k,q,T t )Pi(k ■ q), 



(36) 



3.2 Integral solution 



For given comoving wavenumber and momentum, Eq. ( 26 \ 
is a linear first order ordinary differential equation, and has 
the explicit solution 

*>«/(k,q,Ti) 



<5/(k,q,r) = e 
. m dfo 



g dg 



Vr'Mk, T')dT', (27) 



where t; is some initial conformal time, and we have used 
the variable (sometimes referred to as the "superconformal 
time" ) 

- r dr' 



s(r) 



(28) 



We can now evaluate the average of Sf over the direction of 
momentum q (this will be needed shortly): 

(29) 



Sf(Kl,r) = — I d g<5/(k,q,-r) 



+ m 



Sf (k,g, n,r 
dfo 



dq 



k j T jl ( fc ^( s - s '))a(^')0(k,-r')dr', (30) 



where ji is the order-one spherical Bessel function of the 



first kind (see for example Weisstein 2012b I, and we have 
defined 



Sf '(k,q,Ti,T) = i / <r<i <• 



<5/(k,q,Ti 



(31) 



The mass density of the considered particle can be evaluated 
by integrating the phase-space density over momenta: 



(32) 



Phdm(x, t) =ma Id g/(x, q, r) 



We can now obtain the density perturbation in Fourier 
space, 

op hdm (k, r) = ma~ s J d 3 g<5/(k, q, r) 



ma / 47rg dq <5/(k, g, t) 



(33) 



Dividing by the mean density p hdm (obtained from Eq. ( 32 \ 
with / = fo), using Eq. ( |30[ ) and integrating the second term 
by parts, we obtain the following expression for the density 
contrast <5hd 



>(k,r) 



Phd: 

si 



k 



^hdm (k, Ti , T 

k 



s) 1 



a(r')^(k,T')dr'. (34) 



In the above equation, 5hdm (k, Ti, t) is the value of the den- 
sity contrast evolved from n to r in the absence of any 
gravitational potential, 

x i s _ J dqq 2 Sf I (k,q,Ti,T) 

Ohdm (k, Ti , r ) = =- . (35) 

J d<? <r/o(<?) 



we obtain the following expression for 5 1 : 

Ohdm k,r»,r) = > f , UM ■ (37) 

J dqq 2 f (q) 



The dimensionless function I in Eq. (34 1 is the Fourier trans- 



form of the unperturbed distribution function in momentum 



space, normalised so that 1(0) = 1 ( Brandenberger et al. 
1987[ |Bertschinger fc Watts|1988[ ), 



I[X;f ] 



Jdqj ( q X) q 2 f (q) 
fdqq 2 fo(q) 



(38) 



For neutrinos described by a relativistic Fermi-Dirac dis- 
tribution, we provide an accurate fitting formula for this 
function in Appendix [C| 

Finally, since the gravitational potential is sourced by 
the total matter density through the Poisson equation 



k 2 (f> 



4ttci*p m 5m = ~Ho^-8m, (39) 
where (Jm is the matter fraction at the present day, we can 



rewrite Eq. (34 1 as 

<5hdm(k,r) = <5hdm(k, Ti,r) 
X 



k 

— 0- 

m 



(s-s')5u(Kr')dr'. (40) 



Note that 0h dm at time r depends, in principle, on its value 
at all prior times through 5m(t' < r). 



3.3 Discussion: limiting regimes 

The unperturbed phase-space density is in general charac- 
terised by a typical comoving momentum go, such that /o(g) 
decreases rapidly for g » q . For neutrinos, for example, 
go = T„,o. We also assume that this is the case for Sf; in- 
deed, provided this is true initially, it will remain true as 
long as linear theory is valid, as can be seen from Eq. |27|. 



3.3.1 Large scales 

In the limit kqo(s - Si)/m < 1, we may Taylor-expand the 
spherical Bessel functions and obtain, up to terms of 0(k 2 ), 

Shdm(k ->■ 0,t) = 5hd m (k,Ti) - aj0hdm(k,Ti)(s - Si] 

- k 2 /' T a(r')(s-s')0(k,r , )dT', (41) 



(42) 



where 5tdm(k, Ti) is the initial overdensity: 

, r , Jdqq 2 SMk,q,T l ) 
J dqq 2 fo(q) 

and #hdm(k, r,) is the initial bulk velocity divergence: 

a n \ l i J dqq 2 (q/mai)Sfi(k,q,Ti) 
hdm (k, Ti)=-k — 



(43) 



Jdqq 2 fo(q) 
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Equation (41 1 is nothing but the explicit solution to the 



(Newtonian, sub-horizon) linearised fluid equations in an ex- 
panding universe (see for example |Ma fc Bert schingcr 1995] ): 

9 hdm = 0, (44) 



<3hd 
dhdm 



~ Ohdn 



k 2 cf>. 



(45) 



We therefore recover the well-known fact that on scales much 
larger than their free-streaming length, HDM particles be- 
have as a pressureless ideal fluid. 

In general, the HDM component should behave like a 
cold species on scales much larger than its free-streaming 
length, regardless of whether it is linear. It is therefore im- 
portant, for our linear approximation to work across all 
scales, that the CDM is indeed linear on scales larger than 
the free-streaming length, otherwise we would not get the 
correct long- wavelength behaviour. This is ensured by the 
hierarchy of scales k{ s < k n \. 

3.3.2 Small scales 

Let us define the function 

„2 



£(r,x,q) 



2m 



+ ma 4>{t, x). 



The change of E along a particle's trajectory is 
dE\ dE dx dE . dq dE dE 



_ dE + dx 

dr Itraj dr dr 



dE dq dE 
dx dr <9q 



(46) 



(47) 



where the partial time derivative is at constant x and q and 
we have used the geodesic equations for x, q to cancel out 
the last two terms. We therefore see that 



dE\ 
dr I 



<, aHE, 



(48) 



provided the potential varies on a Hubble timescale. This is 
an upper limit: if ma 2 |<^| <C q 2 /m then the rate of change 
of E can in fact be much smaller than aHE. 

Let us now consider a perturbation with scale small 
enough that the crossing time is much shorter than the 
Hubble time (which is itself of order of or shorter than the 
timescale over which E changes). In that case, the quantity 
E is an adiabatic invariant during the crossing of the pertur- 
bation. We assume that any particle within the perturbation 
can be traced back to a large distance from the perturba- 
tion at some earlier time t;, where ma 2 (f>(Ti,Xi) <C qf/m 
and hence Ei — E ~ q 2 / (2m) (the distance must however 
be small enough to be crossed in a timescale short compared 
to the time for E to change significantly). If furthermore at 
that initial time the phase-space density is nearly unper- 
turbed, /(t;) ~ fo(qi), we deduce, using conservation of 
phase-space, that, for E > 0, 



/o(qu) « MV2m~E) 

fo ( v /g 2 + 2(ma)2<^(r,x)) , 



(49) 



and for E < 0, f(r, x, q) = 0. This argument is generic to 
small-scale perturbations and does not assume linearitjQ 

4 Interestingly, for a relativistic Fermi-Dirac distribution, the 
small-scale overdensity obtained from the linearised version of 
Eq. ( |49| is larger than the one obtained from the full non-linear 
equation. 



If we now assume that \<p\ <C q 2 / (ma) 2 , we can linearise 



Eq. ( 49 1 in (j> and obtain 



xt/ \ ~ {ma) 2 dfo , . 
<5/(T,x,q) « i '——(j)(T,x). 

q dq 



(50) 



We find the resulting overdensity, after integrating Sf by 
parts over momenta, 



where v~ 2 is the mean inverse velocity squared, 



(ma) 



JdqMq) 
Jdqq 2 f (q)' 



(51) 



(52) 



Note, in passing, that these expressions require that 
dfo/dq — > as q — > 0, and would not be valid for a Bose- 
Einstein distribution, for example. 

Therefore, in the regime of shallow potentials, clustering 
on small scales is proportional to the square of the ratio 
of the characteristic velocity dispersion of the potential to 
a characteristic velocity of the HDM particles. Using the 
Poisson equation ( 39 1 , we may also rewrite this in Fourier 
space as 



5hdm(fc) w 2 a 2 H 2 v- 2 n M (a)8M = ( r~ 

2 \K[ a 



S M (k), (53) 



where we have used the definition of the free-streaming scale, 
Eq. (JsJ) and Qm(o.) is again the relative contribution of the 
non-relativistic components to the total energy density at 
scale factor a. This simple dependence is what motivated the 
exact numerical prefactors used in the definition of fef s . This 



result was previously derived by Ringwald & Wong (2004) 



A similar result was obtained for baryon clustering on scales 



much smaller than the Jeans scale by Gnedin & Hui ( 1998 1 



(in that case the Jeans scale plays the role of free-streaming 
scale, even though physically the absence of clustering is due 
to the pressure response rather than free-streaming as is the 
case for neutrinos). 

We could also have used our integral expression Eq. ( |40[ ) 
to reach the same results: if k ^> kf s , then (i) kq/m(s — Si) ^> 
1 and the initial conditions become irrelevant, 8 1 — > 0, and 
(ii) the function X decreases rapidly, such that only the times 
t' rs t matter in the integral. In that case one finds, after a 
change of variables, 



2 r 

> fcfe) w -(ma) (j) I 
Jo 



Xl(X)dX, 



(54) 



which can be shown to give precisely Eq. (51 1 



4 APPLICATION TO SIMULATING MASSIVE 
NEUTRINOS 

Equation ( |40[ ) allows us to compute the density field for hot 
dark matter from a possibly non-linear dark matter poten- 
tial. We shall now specialise to the case where the hot dark 
matter is a neutrino component interacting gravitationally 
with cold dark matter and baryons, whose non-linear evolu- 
tion is computed using an TV-body code. 
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4.1 Initial conditions 

The initial unperturbed distribution function for massive 
neutrinos which decoupled while relativistic is the Fermi- 
Dirac distribution, 



/o(9) = 



9 



1 



ft3 e <,/T„, + l ' 



(55) 



where g — 2, the degeneracy of the neutrino species. Only 
the first two moments of the perturbed distribution func- 
tion at the initial time are relevant: on small scales the ini- 
tial conditions are rapidly forgotten, and on large scales the 
higher-order moments scale as v l and decay rapidly with 
I as neutrinos are non-relativistic. In principle, one should 
extract the full momentum information from a Boltzmann 
code at the startup redshift, and obtain <5/o(k, q,Ti) and 
5fi(k,q,Ti). However, on large scales, neutrinos essentially 
behave as a cold species, with local overdensities and bulk 
flows roughly independent of their individual momenta, and 
on small scales, the exact form of initial conditions is irrele- 
vant. These considerations allow us to assume the following 
simple form for the initial distribution function: 



/„(x, q,Ti) = /o(|q- m I/ a i v„(x, n)\)[l + &,(x, Tj)] 

,dln/ 



fo(q) 



1 + S v (x,Ti) — cnm v v u (x.,Ti) ■ q- 



dq 



(56) 



where we have assumed that the bulk velocity v„ is 
smaller than the characteristic random velocity. Fourier- 
transforming this equation, and using v„(k) = — j;k8 v (k), 
we obtain the multipole moments of our approximate initial 
phase-space density: 



tf/o(k, q,n) = f (q)S v (k,n) 
dfo 



(57) 

8fi(k,q,n) = ^m v k^ 1 a i 9„(k,n), (58) 
5f t (k,q,n) = 0, 1^2. (59) 
The initial conditions propagate to the current redshift ac- 



cording to Eq. (371, which gives us, after integrating the 



second term by parts, 

<5 J (k, n, t) = I Si , s [(5„(k,Ti) - aiO v Qs.,Ti)(a - s,)] , (60) 
where Za 1>aa = Z([s2—si]k/m) and I was defined in Eq. (38 1. 



free-streaming and non-linear scales kf s (z) < k n i(z) for cur- 
rently allowed neutrino masses. For linear scales k < k n \(z), 
the CDM phases are approximately constant and equal 
to the phase at the current time r for any r' < r. For 
scales smaller than the free-streaming scale, the integral in 
Eq. ( |40[ | is dominated by the most recent times, such that 
t' — t\ <S t. It is reasonable to assume that phase evolution 
and mixing are not too important during a short period of 
time, and therefore the phases at all relevant times are ap- 
proximately equal to the phase at the present time. Because 
of the hierarchy of scales k{ s (z) < k n i(z), we can therefore 
assume, for all scales, that the phase and relative amplitude 
at t' < t is equal to that at the current time r, and therefore 



<5 M (k, t') 



Pu(k,r) 



1/2 



<$M(k, 



(61) 



where Pm is the matter power spectrum. We also assume 
that the initial condition piece satisfies a similar relation 



5„(k,Tj) 



Pu(k,T) 



1/2 



5„(k,-r), 



(62) 



and that the ratio 0„(k, Ti)/8„(k, n) = [0 V / 5 v ]i(k) is a func- 
tion of k only. The former is justified because initial con- 
ditions are only relevant on large (and linear) scales, where 



Eq. (62 l is indeed correct, and the latter because structure 



is fully linear at our starting redshift. 



Substituting Eqs. (61 1 and Eq. (62 1 into Eq. (401, with 



the initial condition piece given by Eq. ( 60 1, we obtain, first, 



that the neutrino component is in phase with the CDM, and 
second, that the neutrino power-spectrum is given by 

Py 2 (k,r) = I Sz , 3 Pl /2 {k,n){l-(s-s x )a t [6 v /5Mk)} 
3, 



+ -Q,mH 



Zs>, s PM 2 (k,T')(s - s')dr', (63) 



where I si ,s 2 =I{[s2 — si]k/m). Because neutrino and CDM 
overdensities are in phase (following from our assumptions), 
we can then recover the full neutrino density field from 
P„ 1/2 (fc, r) and the CDM density field by 



5„(k,r) 



P v (k,r) 
p cdm(fc,r) 



1/2 



#cdm(k, r) 



(64) 



4.2 Phase evolution 

Storing the full three-dimensional Fourier transform of the 
gravitational potential at each time step in order to perform 



4.3 Implementation 

We implement the effect of massive neutrinos in Fourier 
space as a modification to the TreePM-SPH code GADGET - 



3 (Springcl 2005 Viel et al. 20101. In order to lower the 



the integral in Eq. (40 1 would require prohibitive amounts of computational cost, and following Viel et al. (20101, we do 



memory. We therefore make some simplifying assumptions 
so that we need only store the power spectrum. 

Previous works ( Viel et al.|2010 Brandbyge & Hannes 



tad 2009 \ used the fully-linear neutrino overdensity (i.e. 



computed assuming even the CDM is linear) and assumed 
that the random phases and relative amplitudes of the neu- 
trino density field did not evolve, but were given by the ini- 
tial conditions. This is true on linear scales, where there is 
no mode mixing, but not on non-linear scales, where neutri- 
nos, even if linear themselves, evolve in the gravitational po- 
tential sourced by non-linear CDM perturbations, for which 
mode-mixing is important. 

Here again, we can take advantage of the hierarchy of 



not compute the short-range tree force due to and experi- 
enced by neutrinos. The size of a PM grid cell for our fiducial 
simulation is 0.5 /i _1 Mpc, corresponding to a wavenumber 
k w 13 h Mpc -1 , a scale two orders of magnitude smaller 
than the neutrino free-streaming length. Neutrino overden- 
sities are therefore completely negligible compared to CDM 
overdensities below the grid scale, and we are justified in 
neglecting the Tree force due to them. However, neglect- 
ing the Tree force experienced by neutrinos is not strictly 
justifiable, and effectively amounts to having a lower resolu- 
tion for the neutrinos than for the CDM ( |Viel et alT]|2010[ >. 
As we shall discuss below, we checked that our results are 
converged with respect to grid size, proving that there is 
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a negligible contribution from scales smaller than our grid 
resolution, where the tree force would be active. 

We evolve CDM particles (and baryons) as well as neu- 
trino overdensities simultaneously as follows: 

(i) We generate adiabatic initial conditions using power 
spectra from CAMB and setting equal random relative ampli- 
tude and phases for the CDM, baryon and neutrino initial 
overdensities. 

(ii) Given the total matter overdensity field up to time 
t — At, we evolve the CDM and baryons with our Tree- 
PM code up to the next timestep r. We then evaluate the 
Fourier-transformed CDM+baryon density, as is required for 
computing the long-range particle forces, and compute the 
CDM+baryon power spectrum at time r. 

(Hi) In order to evaluate the neutrino overdensity at 
time r with Eq. ( |63[ ), we require the total matter power- 
spectrum up to (and including) the current time r. Equa- 
tion ( 63 1 is therefore strictly speaking an implicit equation 
for P„(fc,r), and should be solved iteratively. Since the in- 
tegrand vanishes at t' = t (because of the (s — s') term), 
and neutrinos contribute a small fraction of the dark mat- 
ter, with a good initial guess a single iteration is sufficient. 
We show this in detail in Appendix[B] We therefore approx- 
imate P v (k,r) w P v (k,r — At) to evaluate P v (k, t') inside 
the integral of Eq. ( |63[ ) - we do so by interpolating over our 
tables of stored power spectra. We checked that P„(fc, t) ob- 
tained from this first iteration is converged at the level of 
~ 10~ 4 by comparing with the output of a second iteration, 
which is therefore not required. 

(iv) We then recover the full phase information for the 
neutrino overdensity from Eq. | |64[ ), and evaluate the total 
matter overdensity at the current timestep r from 



(5m (k, T) = (1 - /„)<5cdm,b(k, t) + /„<5 v (k,T). 



(65) 



We finally compute and store the total matter power spec- 
trum at time r, and reiterate steps (ii) to (iv) until the end 
of the simulation. 



4.4 Simulation Details 

The parameters of our simulations are shown in T able[T] O ur 
cosmological parameters were similar to those in |Bird et al.| 
| 2012[ ). We set h = 0.7, the total matter fraction fi M = 0.3, 
the scalar spectral index n s = 1 and the amplitude of the 
primordial power spectrum, A s = 2.43 x 10 -9 . Our initial 



conditions were generated, from transfer functions generated 
by CAMB, using our own version of N-GenlCs modified to use 
2LPT ( |Scoccimarro|1998 \ and account correctly for baryons 
in the initial transfer functiorj^] Where it was important, we 
used fl b = 0.05. 



In a previous work (Bird et al. 20121 we found that a 



significant limitation to the inclusion of massive neutrinos 
was the restriction on initial redshift due to the slightly rel- 
ativistic behaviour of the neutrinos at early times. Although 
our method for computing the growth of neutrino perturba- 
tions is non-relativistic, we can at least easily include the 
correct contribution of neutrinos to the background evolu- 
tion by computing 



Name 


m„ (eV) 


Box (Mpc/i" 1 ) 


CDM 


Zi 


K /3 

oar 


SOO 





256 


512 


49 




S00P 





256 


1024 


49 




S05* 


0.05 


256 


512 


49 




S10* 


0.1 


256 


512 


49 




S15* 


0.15 


256 


512 


49 




S20* 


0.2 


256 


512 


49 




S40* 


0.4 


256 


512 


49 




F10BA* 


0.1 


60 


512 


49 


512 


S05IC 


0.05 


256 


512 


99 




S10P* 


0.1 


256 


1024 


49 




S03NH 


0.03, NH 


256 


512 


49 




S03IH 


0.03, IH 


256 


512 


49 





Table 1. Summary of simulation parameters. m„ is the mass of 
one neutrino species. Rows marked with a * were run using both 
Fourier and particle neutrinos. Cosmological parameters were the 
same for all simulations and are given in the text. Where neu- 
trino particles were included, the same particle number was used 
as for the dark matter particles. Simulations with m„ = in- 
cluded massless neutrinos. Note we held SIm constant, making 
f^cdm dependent on fi„. S03NH and S03IB had the same total 
neutrino mass as S05, but had three species with masses follow- 
ing the normal and inverted hierarchies, respectively. 



PM=a 3 j t{<l,a)fo(q)4:Kq 2 dq, 



(66) 



5 Our initial conditions generator is freely available at |http:/7| 
github . com/sbird/S-GenIC 



where e(q, a) = sj ml + a~ 2 q 2 is the total energy of a neu- 
trino with comoving momentum q. Thus, p„(a) interpolates 
between the behaviours of matter and radiation. This is ex- 
tremely easy to implement with our Fourier method, but 
somewhat harder in a purely particle simulation, where the 
energy density is by default computed from the particles' 
rest-mass and does not account for their kinetic energy. For 
consistency, we also included the effect of radiation density 
in the background evolution, and, for CDM only simula- 
tions, the effect of massless neutrinos. Our particle- based 
simulations still have p„ = ~p v a~ :i . 



5 CODE VALIDATION 

5.1 Convergence tests 

In this Section, we check the convergence of our method with 
respect to particle number and initial redshift. 

Figure [2] shows the impact of changing the initial red- 
shift from Zi = 49 to Zi = 99 for neutrinos of mass m v = 0.05 
eV. The change at low redshift is of order 0.5% and is con- 
sistent with slightly more small-scale non-linear growth for 
the simulation with the higher initial redshift, as expected 
from a pure cold dark matter simulation. 

Figure [3] shows the effect of changing the particle 
number, comparing S10 to S10P, which has 8 times more 
particles. Changing the number of particles causes the 
sampling of initial structure to change slightly, introduc- 
ing sample variance. To minimise this effect we show 
the ratio of the suppression due to neutrinos for pairs 
of simulations with different particle numbers, that is, 
(Psio/Psoo)/(Psiop/Psoop) - 1. Our results are converged 
at the 1% level for k < 10ft" 1 Mpc, and at the 0.2% level for 
k < lh" 1 Mpc. The effect of increasing particle resolution is 
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l.Or 



rs 0.5- 



•% 0.0 



-0.5 





k/(h Mpc-1) 

Figure 2. Change in the total matter power spectrum between 
S05 (initial redshift z t = 49) and S05IC (initial redshift Zj = 99). 
Positive values correspond to more power in S05. Each line shows 
a different redshift: z = 9 (green dot-dashed), z = 3 (grey dotted), 
z = 1 (blue dashed) and 2 = (black solid). 



10" 2 10" 1 10° 10 1 

k/(h Mpcr 1 ) 

Figure 4. The suppression in the total power spectrum at z = 0. 
Solid lines show a 256Mpch —1 with our fourier-based method 
(blue) and the particle method (red). The dashed line shows the 
prediction of linear theory. The spike at k 7 /i/Mpc (corre- 
sponding to half the grid frequency) is a numerical artefact. 



1.0 



0.5 



0.0 



-0.5 



-1.0 




layed onset of non- linear growth, as discussed in |Bird et al.| 
120121). 



10" 



k/(h Mpc-1) 



10" 



Figure 3. Effect of the number of CDM particles on the power 
suppression by neutrinos. What is plotted is the fractional dif- 
ference between the ratio (Psio/^soo) (both using 512 3 CDM 
particles) and the ratio (fsiop/fsoop) (both using 1024 3 CDM 
particles). Positive values correspond to a larger suppression of 
power in the 1024 3 -particle simulation. Each line shows a differ- 
ent redshift: 2 = 9 (green dot-dashed), 2 = 3 (grey dotted), 2 = 1 
(blue dashed) and 2 = (black solid). 



to increase the overall growth of non-linear structure, and 
as a consequence to enhance the relative suppression due to 
massive neutrinos. 



5.2 Total Power Spectra 

Figure [4] shows the suppression in the total matter power 
spectrum, defined as the ratio between the total matter 
power spectrum including massive neutrinos with M v = 
^2 m v =0.3 eV and the matter power spectrum with mass- 
less neutrinos, but unchanged Q.yi. The particle and Fourier 
neutrino simulation methods agree extremely well, and we 
see again the enhancement in the suppression due to the de- 



Figure [5] shows the percentage difference in the total 
power spectrum between the particle and semi-linear neu- 
trino simulations, for redshifts between z — 3 and z — 
and for a variety of neutrino masses. The agreement is ex- 
tremely good in all cases. We show those scales where our 
simulation is resolved to less than 1%; note that the differ- 
ence between the two methods is generally smaller than the 
convergence error shown in Figures [2] and [3] The effect of 
non-linear growth in the neutrino component, neglected in 
our semi-linear method, only begins to become important 
for M„ = 1.2 eV, at z < 1. Since neutrinos of this mass are 
already ruled out by current data, this is unlikely to be a 
practical limitation. 

One interesting feature is that on small scales we seem 
to see (slightly) more power in the Fourier method than 
the particle method, and the effect increases slightly at high 
redshift. This is not due to shot noise, which would lead to 
increased power in the particle method. This appears to be 
a convergence effect; as we have shown in Section [5.1| these 
scales are not converged to less than 1%, and increasing the 
resolution of the simulation leads to a marginal increase of 
power on small scales. We therefore suspect that our semi- 
linear method has an effective resolution slightly higher than 
the particle method, although if we computed the Tree force 
for the neutrino particles, the situation would probably be 
reversed. 

It is interesting to compare our results to those obtained 
using the original fully-linear Fourier-space neutrino method 
of |Brandbyge et all ( |2T)08| ) . Fi gure 1 of that work is compara- 
ble to our Figure J5j Our semi-linear method agrees with the 
particle-based method slightly better in the linear regime. 
This is due to our derivation of the neutrino power spectrum 
from the CDM, which incorporates sample variance in the 
neutrino component in a way very similar to the particle 
method. In the non-linear regime, the agreement between 
the semi-linear method and the particle method is signifi- 
cantly better than that between the fully-linear and particle 
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Figure 5. Percentage differences in the total matter power spectrum between the semi-linear Fourier method and particles. AP t = 
^Fourier — ^Particle: so positive values correspond to more power in the Fourier method. The simulations used were S05 (Top left), S10 
(Top right), S20 (Bottom left) and S40 (Bottom right). Each line shows a different redshift: z = 9 (green dot-dashed), 2 = 3 (grey 
dotted), 2 = 1 (blue dashed) and 2 = (black solid). 



methods. The fully-linear method shows differences from the 
particle method of > 1% for neutrinos with M v = 0.6 eV, 
and has an error of about 5% at 1.2 eV. We ran test simula- 
tions using the fully-linear theory power spectrum, but with 
our improved model for the neutrino phase structure (i.e. 
assuming the neutrinos are always in phase with the CDM 
rather than keeping their initial phases). The result was still 
in good agreement with Brandby ge et ah] ( |2008[ ), showing 
that the improved behaviour of the semi-linear method in 
the non-linear regime does indeed result from accounting 
for the deeper non-linear potential wells sourcing the linear 
growth of neutrinos. 



Bird et al. (20121, using the linear theory neutrino im- 
plementation of Viel et al. (20101, found somewhat larger 



differences between the linear and particle methods than 
|Brandbyge et al.| ( [2008] ). When preparing this work, we 
discovered that this was due to a slight inconsistency in 
the algorithm used; the Fourier space neutrino implemen- 
tation was altering the initial CDM transfer function to in- 
clude neutrinos, adding the neutrino power when calculating 
long-range forces, but neglecting them when outputting the 
power spectrum. Once this was corrected, we obtain fully- 
linear Fourier-space neutrino results in good agreement with 



5.3 Neutrino Power Spectra 

Figure [6] shows the power spectra for the neutrino compo- 
nent at two different redshifts for the different methods. In 
the linear regime both codes agree well with linear theory. 
On non-linear scales our semi-linear method shows addi- 
tional power over the fully-linear neutrino power spectrum. 
On all scales, Pcdm > P v and k 3 P v /(2ir 2 ) < 1. 

At z = 0, the particle method produces additional 
power in the neutrino component. It begins to deviate from 
our semi- linear Fourier based method at k = fc n i, the non- 
linear scale for the dark matter. This first becomes apparent 
for z < 0.5, for all M v ^ 0.1 eV. We have checked that this 
area of the power spectrum is unchanged with a 512 3 par- 
ticle simulation, verifying that it is not being affected by 
shot noise. We ascribe this effect to non-linear growth in 
the neutrinos around the centres of clusters, as described in 
iRingwald & WonglpO 



Villaescusa-Navarro et al. (2011 1, 



those of Brandbyge et al. ( 2008 1 



and as we discussed in Section [2T4] 

This does not appear to have any effect on the total 
matter power spectrum. Most of the effect of neutrinos on 
the total matter power spectrum comes from their time in- 
tegrated evolution; in particular the size of the trough is 
due to their effect on the onset of non-linear growth. Since 
the highly non-linear effects on the neutrino power spec- 
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Figure 6. The power spectrum of the neutrino component at (Left) z = (Right) 2 = 1. Solid black shows the results of our semi-linear 
fourier-based method from simulation S10. Dashed red shows the particle method, from S10P, to obtain lower shot noise. Dotted green 
shows pure linear theory. The vertical dashed grey line shows the approximate non-linear scale for the dark matter. 



trum only take place at relatively late times, the impact on 
the total power spectrum is minimal. Moreover, on scales 
much smaller than the free-streaming scale, the main effect 
of including neutrinos is essentially to reduce the matter 
overdensity by a factor (1 — f v ), while keeping the same 



background expansion rate ( Lesgourgues & Pastor 2006 1 . 
As long as \S V \ <C |<5 c dm|, the exact value of the neutrino 
overdensity does not matter very much. The same conclu- 
sion was reached by |Shoji fc Komatsu ( 2009 1 , who com- 
pared their thrid-order perturbation theory calculations (for 
both CDM and neutrinos, the latter approximated as a per- 
fect fluid with pressure) to a computation where CDM is 
treated to third order but neutrinos are computed to linear 
order (without accounting for non-linear growth of poten- 
tials), similar to the treatment of Saito et al. (20081. They 
found that although neutrino overdensities can be underes- 
timated by order unity on non-linear scales, the total matter 
power spectrum is very accurately obtained, even with the 
simpler methocj^] 

Finally, let us point out that for z > 0.5, subtracting 
the scale-free shot noise from P v in the particle simulation 
produces very good agreement with the semi-linear Fourier 
method, even at scales where the shot noise dominates. This 
is further evidence that neutrino shot noise is not having a 
strong dynamical effect, and is not causing spurious cluster- 
ing. 



5.4 Performance 

Simulations S05-S20 were consistently faster when using our 
Fourier method. The speed increase was 13% of the total 
walltime (which includes time spent reading and writing to 
disc). Note that the slowest single algorithm in GADGET is 
the Tree method for computing short-range forces, which 



B In the notation of |Shoji fc Komatsu| 1 2009 ), our method as- 
sumes Ptot = fc Pqo, T+ (2/ e /„gi + fi g()Poo,c, which is better 
than the treatment of Saito ct al. (2008 ), in the sense that we use 
the full non-linear CDM power spectrum as a source for neutrino 
overdensities. 



is disabled for neutrinos even for our particle based sim- 
ulations, hence a large proportion of the execution time is 
independent of the method used to simulate neutrinos. More 
importantly, the total memory usage of GADGET was 40% 
smaller in the Fourier method than with particle neutrinos, 
essentially identical to the memory usage of a pure dark mat- 
ter simulation. This is important because memory is often 
the limiting factor when performing large modern simula- 
tions. 

The S10P simulation, which had 8 times more dark mat- 
ter particles than S10, took 12 times longer. This scaling is 
similar to that expected for a pure dark matter simulation, 
demonstrating that our neutrino method scales well. In fact, 
the only limit to scalability in the neutrino calculation is the 
need for inter-process communication when computing the 
power spectrum. 

Overall, our Fourier method appears to have similar 
performance characteristics to a pure dark matter simu- 
lation, as should be expected; the time to compute the 
neutrino power spectrum is completely negligible compared 
to the JV-body algorithms, and the most costly part of 
our Fourier algorithm is summing modes on the Fourier- 
transformed density grid to compute the power spectrum. 



6 APPLICATIONS 
6.1 Lyman-Q Forest 

The Lyman-a forest is an indirect probe of the matter power 
spectrum at small, non-linear scales (A; = 0.1 — 4/i -1 Mpc), 
and at high redshift, z = 2 — 4. The power spectrum of the 
Lyman-a flux measures the clustering of the absorption sig- 
nal from neutral hydrogen in quasar spectra, and can be used 
to place constraints on the amplitude of primordial pertur- 
bations. When combined with constraints from large scales, 
this can lead to tight constraints on neutrino mass | |Seljak| 



2000| |Gratton et~aT1|2008| |Viel et al.||2010[ ). At these high 



redshifts, shot noise could be an issue for light neutrinos, so 
it is a natural place to apply our method. In addition, sim- 
ulations with neutrino particles, dark matter and baryons 
can become unwieldy. 
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Figure 7. The change in the flux power spectrum between 
our Fourier method and the particle method, using simulations 
F10BA with ra v = 0.1 eV. Positive values correspond to more 
power with the Fourier method. Each line shows a different red- 
shift z = 4 (blue dotted), 2 = 3 (green dashed), and 2 = 2 (black 
solid). SDSS Lym an-a; data constrains th e flux power spectrum 
for k ^ 4ft" 1 Mpc |McDonald et al.| ( |2006| . 



We used simulation F10BA to simulate the Lyman- 
a forest at z = 2 — 4. Since the Lyman-a forest is not the 
main focus of our paper, we shall not explain in detail our 
simulation methodology, and refer the interested reader to 



Viel et al. (20101. Our simulation calculates the clustering 
and ionisation state of the gas at z = 2 using smooth par- 
ticle hydrodynamics, together with radiative cooling and a 
reaction network which assumes optically thin gas in ionisa- 
tion equilibrium. We then simulate the flux in a set of 16000 
quasar spectra by calculating the absorption along random 
lines of sight. The averaged power spectrum of this flux is 
the observable quantity. Although constraints are available 
for k £ 4/i _1 Mpc, the Lyman-a forest is also sensitive to the 
collapse scale of hydrogen absorbers k ~ 60h~ Mpc, so very 
high resolution simulations are required. Figure [7] shows the 
change in the flux power spectrum between our two meth- 
ods; clearly the differences are quite small. For simulations 
with the resolution of F10BA and m„ = 0.1 eV, therefore, 
shot noise is not having a significant affect on the dynamics, 



even at z = 4, as also found by Viel et al. (2010). The total 



effect of neutrinos on the flux power is 5 — 10%. 



6.2 Hierarchy 

Achieving maximal accuracy for small neutrino masses re- 
quires us to account for the neutrino hierarchy. This is diffi- 
cult in particle-based codes, as two neutrino species doubles 



the amount of memory required, but has been done (Wag- 
|ner et aL] [2bl2). It is, however, essentially trivial using our 
method. We have performed two simulations (S03IH and 
S03NH) with a total neutrino mass of M v =0.1 eV using 
either the normal or inverted hierarchies, both to demon- 
strate the technique and examine the effects of the hier- 
archy. For completeness, the masses of the three neutrino 
species were 0.022, 0.024, 0.054 eV for S03NH and 0, 0.049, 
0.051 eV for S03IH. Our results are shown in Figure [8] to- 
gether with the prediction of linear theory. The linear theory 



Figure 8. The difference in the total matter power spectrum 
between the inverted and normal neutrino hierarchies. Negative 
values correspond to more power in the normal scenario. Each 
line shows a different redshift: 2 = 3 (grey dotted), 2 = 1 (blue 
dashed) and 2 = (black solid). The orange dot-dashed line shows 
the linear prediction at z = 0. 



change due to the hierarchy is quite small (below our con- 
vergence error). Furthermore, the effect in linear theory is 
mostly due to the change in background evolution caused by 
one of the possible hierarchies having one massless species. 
As shown in Figure 16 of Lcsgour gues fc Pastor| ( [2006] ) , the 
difference between the hierarchies when all three neutrino 
species are massive in both cases is much less. Although 
differences in our cosmological parameters mean that the 
linear effect is smaller than that found by | Wagner et al.| 
(2012), we still find, as they did, an enhancement at non- 
linear scales, analogous to the enhancement for the total 
neutrino effect shown in Figure [4] Determining the hierar- 
chy separately from the total neutrino mass from an effect 
this subtle is likely to prove challenging. However, obtaining 
robust constraints with M v ^ O.leV will require accounting 
for the mass splitting. 



6.3 21cm Forest 

Observations of the spin-flip transition of neutral hydrogen 
in the 21cm forest have the potential to probe structure 
growth at z ~ 20 ( |McQuinn et al.|2006[ ) . Structure formation 
at these redshifts is sensitive to the non-linear growth of the 
first, small halos. No-one has yet examined the impact of 
neutrinos on the 21cm forest, and performing the requisite 
simulations is beyond the scope of this paper. However, our 
code would be ideal for such a study. Halos at these redshifts 
are not yet sufficiently massive for our method to break down 
in their central region, so one could probe arbitrarily small 
scales. Furthermore, shot noise would be a severe issue for 
the particle method at these redshifts. 



7 DISCUSSION AND CONCLUSIONS 

We have presented a new semi-linear method for simulating 
the effects of massive neutrinos on cosmic structure, tak- 
ing advantage of the hierarchy between the neutrino free- 
streaming scale and the non-linear scale, k{ s < k n i for cur- 



© 2012 RAS, MNRAS Q00,|ilfT5l 



14 Y. Ali-Haimoud and S. Bird 



rently allowed neutrino masses. The power spectrum of the 
neutrino component is calculated using perturbation theory, 
and added to the cold dark matter in Fourier-space, as in 
Brandbyge et al. (20081. However, we improve upon their 



method by sourcing neutrino perturbations using the grav- 
itational potential obtained from the full non-linear dark 
matter overdensity rather than from the dark matter power 
predicted by linear theory. We evolve the CDM, baryons and 
neutrinos simultaneously and self-consistently. 

For observationally relevant neutrino masses, our 
method gives results for the non-linear power spectrum 
essentially identical to a fully converged neutrino particle 
treatment, a significant increase in accuracy over fully lin- 
ear Fourier-space neutrinos. Our method also has several 
advantages over a particle treatment; it is faster, uses signif- 
icantly less memory, and does not suffer from shot noise in 
the neutrino component. Furthermore, it allows for an easy 
inclusion of the correct background evolution with relativis- 
tic corrections, and makes inclusion of the neutrino hierar- 
chy trivial. These properties make it especially suitable for 
simulating neutrinos at the low end of the currently allowed 
mass range. 

Our treatment is not strictly valid in the inner regions 
of very massive halos, where the escape velocity might be 
larger than the thermal velocity of neutrinos, which may 
be captured and cluster non-linearly. Accurately describing 
the distribution of neutrinos in these halos requires treating 
them as iV-body particles. Because of this effect, we do not 
accurately recover the non-linear neutrino power spectrum 
at redshifts z < 0.5. However, this does not affect the total 
matter power spectrum, because on the one hand most of 
the neutrino effect is imprinted at earlier times, where our 
method is extremely accurate, and on the other hand, even 
with an enhanced clustering due to non-linear growth, neu- 
trino overdensities are very subdominant compared to those 
of the cold dark matter. 

While we have focussed on neutrinos, our treatment can 
also be easily applied to any hot dark matter particle whose 
characteristic velocity is larger than (or at least of the order 
of) the escape velocity of massive clusters, and whose free- 
streaming length today is larger than the non-linear scale. 
There are many applications of our method beyond the non- 
linear matter power spectrum; we have briefly discussed the 
Lyman-Qf and 21cm forests. It allows the inclusion of neutri- 
nos in an iV-body simulation at minimal cost and effort, and 
should thus be useful for any problem in non-linear structure 
formation. 
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APPENDIX A: ANALYTIC APPROXIMATION 
FOR LINEAR HDM CLUSTERING 

In this appendix we derive an approximate expression for 
Shdm valid on all scales during matter domination. We do 
not use this expression when evaluating neutrino clustering, 
but it provides some insight. Our main simplification here is 
to assume that the gravitational potential is nearly constant 
in time (or equivalently that the matter overdensity grows 
roughly linearly with the scale factor). During matter domi- 
nation and for linear evolution, (f> is in fact strictly constant. 
With this approximation, we get 



i(k,r) w im(k, 



(k,r) 



rs — Si 


' k : 


/ 1 


— s 


'0 


m 



[a(s — s)] 2 s ds, 



(Al) 



where we made the change of variables s 

T = 



a(s — s) = 



matter domination, we have a oc r 2 

ajs) 

,1 + ±a 2 (s)H(s)s) 2 ' 

Changing variables to X — (k/m)s, we obtain 

5 h dm(k,r) w <5hdm(k, Ti,r) 



= s — s'. Assuming 
2/{aH), and 



(A2) 



(ma) 0(k, r) 



X1(X) 

+ x/x k y 



rdX, 



(A3) 



where X k = 2k/(ma 2 H) ~ q ( 7 1 A:/fcf s (a). Let us now consider 
scales small enough that initial conditions are "forgotten", 
kqo(s — Si)/m 3> 1. In this case 5 1 « and the upper limit 
of the integral in the above equation can be approximated 
as infinity: 



3hdm 



(k. 



-(ma) 0(k, t) 



XT(X) 



l0 (l + X/X k ) 4 
Note that during matter domination, 
faikqo{s — Si) 



qoX k 



rn 



dX. (A4) 



(A5) 



so provided a/ai is sufficiently large, qoXk may be of order 
unity, even if kqo(s — Si)/m 2> 1. Using Poisson's equation 
( 39 1, we arrive at 



<5 hd m(k, r) « J r (fe//cf s )5 M (k, r), 

where we defined the dimensionless function 



X l J< 



XI{X) 

(i + x/x k y 



T dX. 



(A6) 
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For k < fc fs , q X k -> and T(k < fe fs ) 
scales, <$hdm ~ 8 M - For k » qoX k - 
show that 



_6_ 
Xf. 



F{k > fc fs ) = ^2 q- 



where 



■> 1, i.e. on large 
co and one can 



(A8) 



(A9) 



^2 _ fdqfojq) 
Jdqq 2 f (q) 

In th is limit we recover the small-scale solution of Section 
m = -v~ 2 <f> = (k is /k) 2 5m- 



3.3.2 



We evaluated the function T numerically for massive 
neutrinos and found that it is approximated to better than 
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12% for any value ol k/k& by the very simple form (see also 
|Wong| | |2008l l) 



(A10) 



l(x) 



APPENDIX B: ACCURACY OF THE 
ITERATIVE SOLUTION FOR THE NEUTRINO 
POWER SPECTRUM 

In this appendix we discuss the accuracy of the step where 
we obtain the neutrino power spectrum at time r given the 
total power spectrum at earlier times Pm(j' ^ t — At), and 
the current CDM power spectrum, P c dm(T)- We start by 
noticing that, because of our assumption of total correlation 



of neutrinos and CDM [Eq. |64l], we have 
PU 2 (k,r) = (1 - U)PUl{KT) + f v Pl'\k,r). 



(Bl) 



We now split Eq. 1 63 1 in a term that is completely known 



at the current time step and a term that depends on the 
neutrino power spectrum between r — At and r 



P„ 1/2 (fc, r) = Pj /2 (fc, r; known) + AP„ 1/2 (implicit) 
where 

AP„ 1/2 (implicit) 



(B2) 



-D.MHo.fy 



(s-s')l a ,, s Pi /2 (T')dr'.(B3) 



We see that this renders the equation for P v {k,r) implicit. 
Now, provided HAt = clHAt <C 1, which is obviously the 
case in a iV-body code, the neutrino power spectrum in 
the integral is approximately P v {t') = P„(r)[l + O(HAt)]. 
Switching to the integration variable s' , and using dr' — 
a(r')ds' « a(r)ds' , we obtain, up to corrections of relative 
order 0{HAt), 



API 1 " 2 (implicit) = e (r)P„ 1/2 (fc,T), 
where we have defined 



e(r) 



3 



(s - s')l s , tS ds', 
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where As ~ At /a ~ At/a . Now 1^1 for any value of its 
argument, and therefore 



^n M H^a- 3 j„{Atf 



-J v Q. m {a){HAt) 2 



< f u {HAtf. 
Therefore, whereas the exact solution is formally 



,1/2 



(k,r) 



,1/2 



l + e(k,T)' 



(k, t; known), 



(B6) 



(B7) 



because e = O (/„(//At) 2 ) is small, an iterative solution 
will converge very rapidly, the error after n iterations being 
of the order of e n [where the n = iteration is P v = and 
the (n + l)-th iteration is obtained by using P„ (implicit) 
obtained from the n-th iteration]. Instead of initialising P„ 
with 0, we moreover initialise it assuming P v (t) ~ P v (t — 
At) in the integral. The overall error is therefore extremely 
small, of the order of O (f v (HAt) 3 ) . 




Figure CI. Function T(x) that determines the damping of per- 
turbations due to free-streaming of neutrinos [Eq. | |38| l with 
x = XT vfi ] 



APPENDIX C: FITTING FUNCTION FOR X{X) 

The function X(X) only depends on the dimensionless prod- 
uct x = XT Vt o (we remind the reader that X was defined as 
an inverse comoving momentum). For neutrinos described 
with the relativistic Fermi-Dirac distribution, it can ob- 



tained from the following series expansion ( Bertschingcr & 
Watts|1988| >: 



l(x) 



3C(3) 



(n 2 + x 2 



(CI) 



We also find that T(x) is approximated to very high accuracy 
by the fit 



I&t{x) 



1 + 0.0168 x 1 + 0.0407 i 4 



1 + 2.1734 x 2 + 1.6787 a; 4 - 1811 + 0.1467 x s 
This fit gives the correct asymptotic behaviours 



(C2) 



Z ( a; ) ~ 1 ~d^ x - 1 - 2.1566 x 2 , x<l, (C3) 
2Q (6) 

m x 1 0.2773 , 

1{X) " W^-. »>1. (04) 

The relative accuracy of the fit is \AI/I\ ^ 1% for $J x 2 
and \AT/T\ ^ 3% for ^ x ^ +<x; the absolute accuracy 
is |AX| & 0.07% for < x < +oo. We show the function 1 
in Fig. |C1| We have checked explicitly that modifying our 
code to use the full function instead of the fitting formula 
does not affect our results. 
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